/**************************************************************************/
/*   treba - probabilistic FSM and HMM training and decoding              */
/*   Copyright © 2012 Mans Hulden                                         */

/*   This file is part of treba.                                          */

/*   Treba is free software: you can redistribute it and/or modify        */
/*   it under the terms of the GNU General Public License version 2 as    */
/*   published by the Free Software Foundation.                           */

/*   Treba is distributed in the hope that it will be useful,             */
/*   but WITHOUT ANY WARRANTY; without even the implied warranty of       */
/*   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the        */
/*   GNU General Public License for more details.                         */

/*   You should have received a copy of the GNU General Public License    */
/*   along with treba.  If not, see <http://www.gnu.org/licenses/>.       */
/**************************************************************************/

/* Numerical algorithms, including fast log2 approximations */

/******************************************************************************/
/* Minimax approximation of log2(exp2(x)+1) through 4th order polynomial.     */
/* Coefficients are chosen by intervals: [-1,0],[-2,-1],[-3,-2],...[-61,-60]  */
/* Calculated with the Remez algorithm (max. err. ~2.95e-7)                   */
/* Array order: x^0, x^1, x^2, x^3, x^4                                       */
/******************************************************************************/

#include <assert.h>
#include <math.h>

inline PROB log1plus_minimax(PROB x) {

    static const double mm[61][5] = {
	//{1.0000000000000000000,0.50001046104880397131,0.086736084604588520792,0.00024467657550238233989,0.0015184462591259856854},
		{1.000000294523528266023080641747018118200,5.000137983039909745599685077166821586592e-1,8.674565581413996689083361220979003743259e-2,2.535058752268224328526728558807895816915e-4,-1.516439960822520490434771183446816832550e-3},
	{1.001540902819038643868034951278937150423,5.053932724630792281026719430573953839270e-1,9.398019953671193331596899635167916987335e-2,4.736317746519460505191161205702483284235e-3,-4.285606238356233021835429915032331912140e-4},
	{1.013976378001803373719080130261091407664,5.302196375839637865890842150875544468037e-1,1.128036571585725980940935215767100213248e-1,1.116060647500509192610719471562426266478e-2,4.038372351507354482020017201054796713189e-4},
	{1.023131887741515939659513867015413784972,5.436431390241952723476902358142572813042e-1,1.201726082355714224040703285933379155073e-1,1.295614207390065258466470987602678285170e-2,5.677144309040374580981548939465879007012e-4},
	{9.835731010812080088792718903383250177176e-1,5.050412800066874725441253191012299753141e-1,1.060138603436639294902261920551732613200e-1,1.064233640714853460177369235785694927148e-2,4.255568505873053481063633077086189126616e-4},
	{8.785837755469130268379738337307369595988e-1,4.212492133056919924809435996347298762156e-1,8.088476028538132830662281420466876267366e-2,7.286063010441202498156555536569382213168e-3,2.571123564851029852691617774209338405663e-4},
	{7.308027596809948836305626723755421141469e-1,3.224898735485838972335374974225996862103e-1,5.609931139717292602758423280134301501202e-2,4.517479281451179182946683596972370600950e-3,1.409757519050094607255642771160615181423e-4},
	{5.729251353692974726328883588237432713562e-1,2.319262332531495481320082839569372770976e-1,3.659728393156119566896355524304121167569e-2,2.649031393921563102265002971696813570262e-3,7.377617211209817543724534818125932575361e-5},
	{4.284662174687129745712318434931123081573e-1,1.593916257690278197207690411432013381941e-1,2.292850909254315600082871464616439604850e-2,1.503305730591576403520118746128495161117e-3,3.773419922889817282851732289801761767787e-5},
	{3.085841212758812309999120320771806912415e-1,1.058863859189567602616804347395344212425e-1,1.396770609670425039470453848858386200084e-2,8.358997265215243173591048436061557046093e-4,1.908170542603252420178862540366703783811e-5},
	{2.155264123882983098749855991600247817735e-1,6.851233959460410183724213278971170157251e-2,8.335958422541212102619996563177562514459e-3,4.585396018125127919513712505533770682245e-4,9.594893492857594836122269201147926166984e-6},
	{1.467306245726877666077626942676748333881e-1,4.339920189972891202144505870112213624964e-2,4.896769013239183127472817654894797060186e-3,2.491225574370352032021045834479735311566e-4,4.811005900201946789976230286497940717666e-6},
	{9.774597747642570611623642585503897029192e-2,2.701135701271180361834345645433818988774e-2,2.840074071729002685693910581519191512437e-3,1.343629108405785720060402064864616879042e-4,2.408898876803461995664849180086577699127e-6},
	{6.390240722468234336692898585104944560666e-2,1.656196793982784935350957620872460454039e-2,1.629841988137060886033206442343713715054e-3,7.204759716181586073189334267641694645559e-5,1.205299188513323577176745175008042077785e-6},
	{4.109530322632024551694888123780769985962e-2,1.002432660091042550516497862044328551315e-2,9.269044219326754916373345637105467655215e-4,3.844733846613674964123871743259133727271e-5,6.028621279112513484528526956144178930942e-7},
	{2.604651704796624806221147750492181508942e-2,5.998857299259351027129693004774952666181e-3,5.230158820585189950977248989959196321287e-4,2.043284216808469960702744546826352735798e-5,3.014842093896247195945892201713717039882e-7},
	{1.629574371510135729387971140997985649272e-2,3.553969682255529864547345125136200672723e-3,2.930853449919225606372295393363827228651e-4,1.082030492014435790730850446793061768234e-5,1.507553925561383747207100592928904011798e-7},
	{1.007725466429307457928882218163015068090e-2,2.086682932279709817708596950126145066405e-3,1.632320249339806286795112234767520211309e-4,5.711905400890433438404547803320168557688e-6,7.538101843127161929560375698734849797613e-8},
	{6.166592810852021628304429555157403713613e-3,1.215315805733406602003594276829073890835e-3,9.041185749446847909978920917266412161871e-5,3.006778600170720770240080230296757498600e-6,3.769133977742138665652662583764853031103e-8},
	{3.737717462473156669258106006521450985779e-3,7.026622053283058117984494282796753011080e-4,4.982968078456188644280170676112499222522e-5,1.578788776152384893993287929005185189189e-6,1.884587753209277764891565287338124725057e-8},
	{2.245914094684305725770080749857493614256e-3,4.035686498423152454763172940607840796264e-4,2.733970162835969618313251809073060759148e-5,8.270905499426057195614983327293335232646e-7,9.422990677258875472042654296680055144308e-9},
	{1.338832664717721788201914319279026374783e-3,2.303840798029541781836884698414285944746e-4,1.493879410404055601672555303086510200499e-5,4.323924099416235579068953347168974922354e-7,4.711508316478430363152009169188885838433e-9},
	{7.922972702672405687699065674711094490131e-4,1.307890086237811189433677560890532741354e-4,8.132130699807709405663734301306179299696e-6,2.256195229860954975265626739800250735402e-7,2.355757402707197225844745208430735781749e-9},
	{4.657234731346026577315424097000851831947e-4,7.386982204094593098754087091647753557274e-5,4.411564764478431250529992566228189279488e-6,1.175213548885805978181224909792554544597e-7,1.177879512471310860977304333925860708915e-9},
	{2.720618620538103096063870592449453738293e-4,4.152512662250294235386843582938633889200e-5,2.385598827740929413533023091255317898203e-6,6.111645692788919966703913409146894898543e-8,5.889399590151730846916027455075052333187e-10},
	{1.580171705274963741730098143850463132205e-4,2.324101836871397046909558051494154234401e-5,1.286241128445408128431121949650261722504e-6,3.173611369941969735477883338518226010925e-8,2.944700302024771549243005586275298086280e-10},
	{9.128823729322136600799461699205922389946e-5,1.295494444820677490156254807866746738526e-5,6.916082012530890839266765652063607943806e-7,1.645699829017103252966042207117727300264e-8,1.472350277749626304855021019882936570645e-10},
	{5.247569910590142135046793616445142420633e-5,7.194060678289050034692570903816580651762e-6,3.709313182968415697803559533749415210589e-7,8.522969558323242550855681945208180816992e-9,7.361751705591250348344593703813574327469e-11},
	{3.002464442449540456321812306316230279927e-5,3.980893426020699241614018914492062186700e-6,1.984709700980208918296314748259221759362e-7,4.408719905863830262694777931537481138819e-9,3.680875932006407067500480958329179247384e-11},
	{1.710422734040868762598444639454100270213e-5,2.195604402388449321176308508239591920320e-6,1.059589922716342243861163654158461109325e-7,2.277977495511747414616719937480764499477e-9,1.840437985805899280496144270537697994597e-11},
	{9.704043605667491263682908071937716380237e-6,1.207214974505590124014617643161619508216e-6,5.645167580837769102123052270029714822467e-8,1.175797513654944728221520532871198877702e-9,9.202189978536236111112806956050839684929e-12},
	{5.484447641224576203381726597065529449632e-6,6.618412653703452350503267709283712470805e-7,3.001714082174804164152957531228620307326e-8,6.063031383797898403328046177453311582210e-10,4.601095001644803025435340889417844059324e-12},
	{3.088458479604903390689693188658224830872e-6,3.618564308589562780175197523528188100416e-7,1.593182842394893138821999205022110470230e-8,3.123537596043724628916392264442200926982e-10,2.300547503916572760527519314184485226062e-12},
	{1.733280697661272233287056762051380371317e-6,1.973331757120187542854428687997321698001e-7,8.441346499330879170202508969316157171779e-9,1.607779749159105324740126938100375879319e-10,1.150273752731829192883736189996345403831e-12},
	{9.696085743590292173102560326724784957429e-7,1.073514018990631532376476703955306105438e-7,4.465291034737605818721360395243583273515e-9,8.268953499066525632385414274164992598604e-11,5.751368765593002996803015177258986392131e-13},
	{5.407542660657723692017343745243259147306e-7,5.826617656970987139962596896850966751768e-8,2.358405230865229758989818177295573682440e-9,4.249504125545619739129738789268168183123e-11,2.875684383279965756601890257363307551726e-13},
	{3.007108152612739136787011420351694731541e-7,3.155581121697607647257436893327796051685e-8,1.243807882731297177491852699085837402269e-9,2.182265750618368043048440349185798101264e-11,1.437842191760848942864078184856755692138e-13},
	{1.667662004078162442507477750365196197228e-7,1.705473504659259956951799930542341586484e-8,6.550692803089486884689982710400498634480e-10,1.119889719190599375617492353842699836805e-11,7.189210959106409875744520300962553590775e-14},
	{9.224363776380102440804314297610964588681e-8,9.199378933794796380053897243067085195251e-9,3.445486622780816285555324388921997055384e-10,5.743232815253642665309746488612302495340e-12,3.594605479628746228230329586681189365969e-14},
	{5.089667226986722137018047548892191372789e-8,4.952924870557156748804625501920938071675e-9,1.809970185281439373218627388776924923170e-10,2.943508517249781067550316354298209946359e-12,1.797302739833258436704936788894355876990e-14},
	{2.801677782038516344958355165843046666879e-8,2.661910662650634497109843058584251325346e-9,9.496769112262315906776942154978907353707e-11,1.507700313429340758317915941323601394079e-12,8.986513699213505489999432228985914748447e-15},
	{1.538758643053035197164792113487652291561e-8,1.428202545949039704456430286944410706895e-9,4.977235557267952106275106755885283243966e-11,7.718231841150908673976030220108117616797e-13,4.493256849618556071618442972575207985175e-15},
	{8.433168824256833527164358422950544309919e-9,7.650403498380251132242583942521669300656e-10,2.605739233309433778235063265261829363000e-11,3.948981057572927184833195256132968137216e-13,2.246628424812228867463908174177381255508e-15},
	{4.612331855583893301257726915641351766482e-9,4.091744076678478527810976190197040857477e-10,1.362778321046620303164709718510336753423e-11,2.019423097284013109685626573804999905097e-13,1.123314212406852141645626380732286770627e-15},
	{2.517668556043846568656183733363678477586e-9,2.185201471374758216019956406730475772415e-10,7.120175012465187609793284679881424765177e-12,1.032177832890477209734800372840581854776e-13,5.616571062036104978012313414859363994793e-16},
	{1.371706328817302179516017326729296854826e-9,1.165361985703673033313902176197454303430e-10,3.716599152485367408168153237766081049177e-12,5.273220585693960621079776522085183661897e-14,2.808285531018513556452202182365282005903e-16},
	{7.460060697873398435554690096701854156388e-10,6.206435840356549075023860180233620003688e-11,1.938240370687552927753139100163725356142e-12,2.692776003467568310554774879654270509258e-14,1.404142765509372045087612472058555689217e-16},
	{4.050178683679676601049722210463806880827e-10,3.301109204107670711766303743756303545759e-11,1.009933068225483015476573048333171694232e-12,1.374470857044027181069215449731212005879e-14,7.020713827547148392591840827660083516925e-17},
	{2.195263541964935173326704258324300343849e-10,1.753623656589639180635522131887936333122e-11,5.257942183832388049473242841197281147524e-13,7.012768561771220717094641254757402974707e-15,3.510356913773646238084365032569057820569e-17},
	{1.187977164264531756608102176993907564716e-10,9.304501861312458327090099781057079168363e-12,2.735215727416921691603014502355565059071e-13,3.576591419161119464701074172078644588424e-15,1.755178456886841129489293671206431650110e-17},
	{6.418966700185583839441433414595898747763e-11,4.931172494095824967894258963597403568756e-12,1.421783288532950779732514753766609368933e-13,1.823399278718305780392810632085263242561e-15,8.775892284434250673564246243633351005494e-18},
	{3.463242499998806440499343947727783680083e-11,2.610517226603860305329831571689912336996e-12,7.385059102157848463640530384474497954671e-14,9.292514239280237429327735045806574130589e-16,4.387946142217136593311567593754028179654e-18},
	{1.865886322849180932233567964249289031758e-11,1.380511857351686805136668728554287345124e-12,3.833233648510797717960088518823766248363e-14,4.734016042484467437568492559634455175235e-16,2.193973071108571110788144914865984070075e-18},
	{1.003909150895296791768415279408779082257e-11,7.293027555134667158634841445482412223363e-13,1.988285256813999683162194355741314271697e-14,2.410887482664406666900931210509709051358e-16,1.096986535554286258927162736930813506477e-18},
	{5.394259651384689949947249921268491380791e-12,3.848980574203442428719130546309148548478e-13,1.030635036607632549333422739749515316323e-14,1.227383472043289447139570270887398625483e-16,5.484932677771433053468539383399344931458e-19},
	{2.894793673005791751968978325760269683645e-12,2.029405835835905259056412377624390019046e-13,5.338928183647988427734040428118065936141e-15,6.246616013771876885631952926327454954548e-17,2.742466338885716966442451116386082280555e-19},
	{1.551567962589901054792803657543910042977e-12,1.069034677089176149199000721561860680068e-13,2.763986071932238298243861608114844112119e-15,3.178157333663653033758089995793957728024e-17,1.371233169442858593148270914364654902117e-19},
	{8.306336675337022760217091368822578009997e-13,5.626366653302989129598511682678622639591e-14,1.430076765921906858110442517105834100357e-15,1.616503330220683752744596997321618529356e-17,6.856165847214293240558967962252323050348e-20},
	{4.441718222138073850369447942491187321448e-13,2.958629470570710761933540111534370365459e-14,7.394916178896801286037976165702941221125e-16,8.219639968047704791288806932600800180116e-18,3.428082923607146688983887328733425427010e-20},
	{2.372529312291007305416992289952648321370e-13,1.554503699235377778851416441800192491369e-14,3.821781113846198396276457417988564726094e-16,4.178381642495995370799785893038057035617e-18,1.714041461803573361668044501268528909826e-20},
	{1.265921709274481217042301748631126620573e-13,8.160998460854488837838924972198941773707e-15,1.974080493999080210366980627204357045868e-16,2.123471650484069163148941014954164434567e-18,8.570207309017866851280474598597185316012e-21}};
    int ptr;
    double xsq;
    ptr = -(int)(x);

    /* Estrin form of polynomial: ex^4  + (dx + c)x^2  + (bx + a) */
    xsq = x*x;
    return(mm[ptr][4]*xsq*xsq + (mm[ptr][3]*x+mm[ptr][2]) * xsq + (mm[ptr][1] * x + mm[ptr][0]));

    /* Horner form (although it has fewer multiplications it is slower on                           */
    /* most architectures due to CPU pipelining issues): a + x (b + x (c + x (d + ex)))             */
    /* return mm[ptr][0] + x * (mm[ptr][1] + x * (mm[ptr][2] + x * (mm[ptr][3] + x * mm[ptr][4]))); */
}

/*******************************************************/
/* Fast log2(1+2^x) table lookup routines for log sums */
/* No error checking: it is up to the calling function */
/* to make sure that 0 >= x >= -LOGTABLESIZE           */
/*******************************************************/

#define LOGPLUSTABLETICKS 128 /* Number of ticks per 1 */
#define LOGTABLESIZE 64

/* Init table, pointer is stored in global var *logplustable */

PROB *g_logplustable;
void log1plus_init(void) {
    int index;
    PROB val;
    long double f;
    g_logplustable = malloc(sizeof(PROB) * LOGTABLESIZE * LOGPLUSTABLETICKS);
    for (index = 0; index < LOGTABLESIZE * LOGPLUSTABLETICKS; index++) {
	f = ((long double)index) / -LOGPLUSTABLETICKS;
	val = (PROB) log2l(1 + exp2l(f));
	g_logplustable[index] = val;
    }
}

/* Straight table lookup, no interpolation */
inline PROB log1plus_table(PROB x) {
    return(g_logplustable[-(int)(x * LOGPLUSTABLETICKS)]);
}

/* Interpolated table lookup */
inline PROB log1plus_table_interp(PROB x) {
    int index;
    PROB val1, val2, w ;
    index = -(int)(x * LOGPLUSTABLETICKS);
    w = -(x * LOGPLUSTABLETICKS) - index;
    val1 = g_logplustable[index++];
    val2 = g_logplustable[index];
    return(val1 + w * (val2-val1));
}

void log1plus_free(void) {
   free(g_logplustable);
}

double *taylor0;
#define TAYLOR_STEP  32
#define TAYLOR_LIMIT 64

void ln1plus_taylor_init(void) {
    /* Precalculate table w/ 2nd order Taylor coefficients for log(2^x+1) */
    /* Uses long doubles for precalculation which are cast into double    */
    int i;
    long double v,n,x2term,x1term,x0term;
    taylor0 = malloc(TAYLOR_LIMIT*TAYLOR_STEP*3*sizeof(double));

    for (i = 0; i < TAYLOR_LIMIT*TAYLOR_STEP; i++) {
	n = -((double)i)/TAYLOR_STEP;
	v = exp2l(2*n) + exp2l(n+1) + 1;
	x2term = 1/(4 + 4 * coshl(n));
	x1term = 1/(1+expl(-n)) - 2*n/(4 + 4 * coshl(n));
	x0term = log1pl(expl(n)) - n/(1+expl(-n)) + n*n/(4 + 4 * coshl(n));
	*(taylor0+i*3+2) = (double) x2term;
	*(taylor0+i*3+1) = (double) x1term;
	*(taylor0+i*3+0) = (double) x0term;
    }
}

double log1plus_taylor(double x) {
    double result;
    int ptr;
    ptr = ((int) -(x*TAYLOR_STEP-0.5)) * 3;
    result = x * (*(taylor0+ptr+2)*x + *(taylor0+ptr+1)) + *(taylor0+ptr);
    return(result);
}

void log1plus_taylor_init(void) {
    /* Precalculate table w/ 2nd order Taylor coefficients for log(2^x+1) */
    /* Uses long doubles for precalculation which are cast into double    */
    int i;
    long double v,n,x2term,x1term,x0term;
    taylor0 = malloc(TAYLOR_LIMIT*TAYLOR_STEP*3*sizeof(double));

    for (i = 0; i < TAYLOR_LIMIT*TAYLOR_STEP; i++) {
	n = -((double)i)/TAYLOR_STEP;
	v = exp2l(2*n) + exp2l(n+1) + 1;
	x2term =  (exp2l(n-1)*M_LN2)/v;
	x1term = -(exp2l(n)*n*M_LN2)/v - (exp2l(n)/v) - (1/v) + 1;
	x0term =  (exp2l(n-1)*n*n*M_LN2)/v + exp2l(n)*n/v + n/v + logl(exp2l(n)+1)/M_LN2 - n;
	*(taylor0+i*3+2) = (double) x2term;
	*(taylor0+i*3+1) = (double) x1term;
	*(taylor0+i*3+0) = (double) x0term;
    }
}

/* Mark Johnson's digamma approx. */
inline double digamma(double x) {
    double result = 0, xx, xx2, xx4;
    if (x == 0) {
	return SMRZERO_LOG;
    }
    for ( ; x < 7; ++x)
	result -= 1/x;
    x -= 1.0/2.0;
    xx = 1.0/x;
    xx2 = xx*xx;
    xx4 = xx2*xx2;
    result += log(x)+(1./24.)*xx2-(7.0/960.0)*xx4+(31.0/8064.0)*xx4*xx2-(127.0/30720.0)*xx4*xx4;
    return result;
}

/* Approximation of inverse chi-square (Weiss and Greenhall, 1996) */
/* Restrictions: df >= 1 and 0.005 <= p <= 0.995                   */
/* df can in theory be nonintegral, but we define it as int here.  */
/* Returns, given the desired p-value and a degree of freedom df,  */
/* the corresponding chi-square value.                             */
/* Max. error of approximation is 3%                               */

double chiinv(double p, int df) {
    double a, A, G, g, y, u, c1, c2, c3, c4, c5, x;
    double p1, a0, a1, b1, b2, t, X, s, b ;
    int k, n, i;
    p = 1-p;
    if (p <= 0.5 && df <= 10) {
	c1 = -0.5748646; c2 =  0.9512363; c3 = -0.6998588; c4 =  0.4245549; c5 = -0.1010678;
	a = (double)df / 2;
	n = (int) a;
	y = a - (double)n;
	G = 1 + y * (c1 + y * (c2 + y * (c3 + y * (c4 + y * c5))));
	for (k = 1; k <= n ; k++) {
	    G = G * (y + (double)k);
	}
	A = p * G;
	u = 0;
	for (i = 0; i <= 7; i++) {
	    g = 1 + (u/(a+1)) * ( 1 + (u/(a+2)) * (1 + (u/(a+3)) ));
	    u = pow((A * exp(u)/g), 1/a);
	}
	x = 2 * u;
    } else {
	a0 = 2.30753; a1 = 0.27601; b1 = 0.99229; b2 = 0.04481;
	p1 = fmin(p, 1 - p);
	t = sqrt(-2 * log(p1));
	X = t - (a0 + a1 * t)/(1 + b1 * t + b2 * t * t);
	s = ( (p-0.5) > 0 ) - ( (p-0.5) < 0 ); /* sign(p-0.5) */
	b = 2/(9 * (double)df);
	x = (double)df * pow((1 - b + s * X * sqrt(b)),3);
    }
    return (x);
}

#include <math.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_sf_gamma.h>
#include <gsl/gsl_cdf.h>

double approx_binomial_test(double p, int successes, int failures) {
    double z, x;
    z = ((double)successes/(double)(successes+failures) - p)/(sqrt(p*(1-p)/(double)(successes+failures)));
    x = 1 - fabs(gsl_cdf_ugaussian_P(z) - gsl_cdf_ugaussian_P(-z));
    return(x);
}

double exact_binomial_test(double p, int successes, int failures) {
    int N, i;
    double p0, p1, logfactsize, logp, logq, mass;
    N = successes + failures;
    logp = log(p);
    logq = log(1-p);
    logfactsize = gsl_sf_lnfact(N);
    p0 = successes * logp + (N-successes) * logq + logfactsize - (gsl_sf_lnfact(successes) + gsl_sf_lnfact(N-successes)); 
    for (i = 0, mass = 0; i < N; i++) {
	p1 = i * logp + (N-i) * logq + logfactsize - (gsl_sf_lnfact(i) + gsl_sf_lnfact(N-i));
	if (p1 <= p0) {
	    mass += exp(p1);
	}
    }
    return(mass);
}

/* Use the test np +- STDDEV * sqrt(npq) which should return      */
/* values between 0 and n for binomial approximation to be valid. */

double flex_binomial_test(double p, int successes, int failures) {
    double n, t1, t2;
    n = (double)successes + failures;
    t1 = (n * p - 5 * sqrt(n * p * (1-p)));
    t2 = (n * p + 5 * sqrt(n * p * (1-p)));
    if (t1 < 0 || t1 > n || t2 < 0 || t2 > n) {
	return (exact_binomial_test(p,successes,failures));
    } 
    return (approx_binomial_test(p,successes,failures));
}

/* Do an MC multinomial goodness-of-fit test on the data in the array c */
/* The corresponding probabilities of H_0 are given in the array p      */
/* Returns: p-value                                                     */

double monte_carlo_multinomial_test(unsigned int *c0, double *p, int size, int iterations) {
    int i, iter, less, more, c0sum;
    unsigned int *c1;
    double p0, p1, *logpi, logfactsize;
    gsl_rng *rng = gsl_rng_alloc(gsl_rng_taus);
    gsl_rng_set(rng, time (NULL) * getpid());
    less = 0; more = 0; c0sum = 0;
    c1 = malloc(size * sizeof(int));
    logpi = malloc(size * sizeof(double));
    p0 = logfactsize = gsl_sf_lnfact(size);
    for (i = 0; i < size; i++) {
	logpi[i] = log(p[i]);
	c0sum += c0[i];
	p0 += c0[i] * logpi[i] - gsl_sf_lnfact(c0[i]);
    }
    /* Generate array c0 from multinomial p      */
    /* iter times, and calculate the probability */
     for (iter = 0; iter < iterations; iter++) {
	gsl_ran_multinomial(rng, size, c0sum, p, c1);
	p1 = logfactsize;
	for (i = 0; i < size; i++) {
	    p1 += c1[i] * logpi[i] - gsl_sf_lnfact(c1[i]);
	}
	if (p1 <= p0) {
	    less++;
	} else {
	    more++;
	}
    }
    free(c1);
    free(logpi);
    return((double)less/(double)(less+more));
}
